!--------------------------------------------------------
! $Id: x2y2model.shar,v 1.2 2011/05/14 21:34:35 daw Exp $
! Description: do GK average of delta_J(0) delta_J(t) vs. t 
! Author:  Murray S. Daw <daw@clemson.edu>
! Created: 14 May 2011
!--------------------------------------------------------
! loop:
!    do some mc steps
!    then do some md steps and calculate delta_J(t) 
! overloop, take average of delta_J(0) delta_J(t)
! delta_J = J - Jbar
! write that out
!
!


program run_gk

  use prec
  use Solve_NonLin
  use md_vv
  use emc
  use global_data
  IMPLICIT NONE

!! controls
!  integer(i12) :: nsamples = 100       ! # of GK samples
!  integer :: nmcsteps = 100       ! # of mcsteps in each sample
!  integer :: nmcequil = 100       ! # of mcsteps to equilibrate (before averaging)
!  real(r8) :: tottime = 1.0        ! how long to run each MD
!  real(r8) :: dt = 1.e-3           ! md time step
!  real(r8) :: snaptime = 1.e-1    ! time between snapshots (J-grid)
!  real(r8) :: du = -1            ! magnitude of attempted displacements (<0, guess)
!  real(r8) :: dv = -1            ! magnitude of attempted changes in momentum (<0, guess)
!  integer :: seed = -123456789    ! ran # seed

  namelist /controls/ lambda, lambda1, nsamples, nmcsteps, nmcequil, &
       &  tottime, dt, snaptime, du, dv, seed

! more vars
! state is defined by 2 vectors (us & vs) which are displacments & momenta
  real(r8), dimension(2) :: us, vs
!  real(r8), dimension(2) :: uscopy, vscopy

  real(r8) :: jtot, j1tot, jtotbeg, jtotend               ! jtot = total curent
!  real(r8), dimension(:), allocatable :: jcorrel, j1correl !holds autocorrelation on each trials
!  real(r8), dimension(:), allocatable :: jauto_avg, j1auto_avg   ! holds autocorrelation average on ten trials

  real(r8) :: ebegin       ! beginning value of energy
!  real(r8) :: emccurrent   ! current value of energy as calculated by emc routines
  real(r8) :: eend         ! beginning value of energy
  real(r8) :: deerr        ! accumlated error

  integer(i12) :: isample         ! sample index
  integer(i12) :: msample         ! sample index
  integer :: nmdsteps
!  integer :: nmdstepsbetweensnaps
!  integer :: nsnaps

  real(r8) :: j0, j1, t0, t
  integer :: isnap, i
  real(r8) :: j02, jjint, taugk

! for action-angle variables
  real(r8) :: action

  real(r8) :: chi
  real(r8) :: num
  real(r8) :: den, den1


!========== Begin Executable Section ==========!

  print *,' <<<*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*>>>'
  print *,' <<<*[]*[]*[]*...                              '
  print *,' <<<*[]*[]*[]*...    run_gk $Revision: 1.2 $   '
  print *,' <<<*[]*[]*[]*...    m. s. daw / 2K8           '
  print *,' <<<*[]*[]*[]*...    daw@clemson.edu           '
  print *,' <<<*[]*[]*[]*...                              '
  print *,' <<<*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*[]*>>>'

! Read parameters from file

  open(unit=10,file='run_gk.controls',status='old')
  read(10,controls)
  close(10)

  write(6,controls)

!
! define initial state (and allocate vars for it)
!
  us = 0.
  vs = 0.

  ebegin = (us(1)**2 + us(2)**2 + vs(1)**2 + vs(2)**2)/2. + &
	      0.25*lambda*(us(1)**2 + us(2)**2)**2 + &
              1./3.*lambda1*(us(1)**3-3*us(1)*us(2)**2)

  print *,' beginning energy = ',ebegin

  nmdsteps = ceiling(tottime/dt)     
  tottime = dt*nmdsteps
  print *,' tottime (possibly adjusted) = ',tottime

  nmdstepsbetweensnaps = ceiling(snaptime/dt)    
  snaptime = dt*nmdstepsbetweensnaps
  print *,' time between snapshots (possibly adjusted) = ',snaptime
  nsnaps = ceiling(tottime/snaptime)
  tottime = snaptime*nsnaps
  print *,' tottime (possibly readjusted) = ',tottime
  print *,' # of md steps between snaps = ',nmdstepsbetweensnaps
  print *,' # of snaps = ',nsnaps

!  print *,'Choose the calculating method: 1--long MD&fft; 2--medium MD&fft; 3--samples_only'

!  read *, i
!  if (i==1) then 
!     call mc_md (us,vs, 1)
!  else if (i==2) then 
!     call mc_md(us,vs,nsamples)
!  else if (i == 3) then 
     call mc(us,vs)
!  else
!     print *,'Invalid input'
!     stop
!  end if


! write a file with correl fn
!  open(file='correl2.dat',status='unknown',unit=10)
!  do isnap=0,nsnaps
!     write(10,9213)isnap*snaptime,jcorrel(isnap)/jcorrel(1)
!9213 format(1x,2g24.16)
!  enddo
!  close(10)

  print *, "lambda = ",lambda
  print *, "lambda1 = ",lambda1

  call closeemc

  stop
end program run_gk


